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We perform a thorough analysis of the relationship between discrete and series representation path 
integral methods, which are the main numerical techniques used in connection with the Feynman- 
Kag formula. First, a new interpretation of the so-called standard discrete path integral methods 
is derived by direct discretization of the Feynman-Kag formula. Second, we consider a particular 
random series technique based upon the Levy-Ciesielski representation of the Brownian bridge and 
analyze its main implementations, namely the primitive, the partial averaging, and the reweighted 
versions. It is shown that the n — 2^ — 1 subsequence of each of these methods can also be 
interpreted as a discrete path integral method with appropriate short-time approximations. We 
therefore establish a direct connection between the discrete and the random series approaches. In the 
end, we give sharp estimates on the rates of convergence of the partial averaging and the reweighted 
Levy-Ciesielski random series approach for sufficiently smooth potentials. The asymptotic rates 
of convergence are found to be 0{l/n^), in agreement with the rates of convergence of the best 
standard discrete path integral techniques. 
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I. INTRODUCTION 

Ever since their introduction over fifty years ago, path 
integral methods have been an intense research field for 
physicists, chemists, and mathematicians alike, even if 
these researchers have sometimes used arguments of a 
rather different nature. The field began with Feyn- 
man's observation that the time propagator of the 
Schrodinger equation can be represented as a "sum over 
histories," effectively giving a formula for the propagator 
as a limit of integrals over spaces of increasing dimension 
In mathematical terms, the existence of this limit is 
problematic though several research directions are known 
k 4, 5, 6]. 

In a significant development, Kag noticed that the 
"imaginary time" version of the formula utilized by Feyn- 
man has a definite probabilistic sense and could be inter- 
preted as an integral of a functional of the seemingly 
ubiquitous Brownian motion 0. Such a formula could 
represent the Green's function for a certain class of diffu- 
sion processes, as for instance the density matrix for the 
Bloch equation. The end product of their work is beauti- 
fully summarized by what is now called the Feynman-Kag 
formula H 



Pfp{x,x';/3) 



Xr{u) + aBl 



du 



(1) 

where p{x, x' \ /?) is the density matrix for a monodimen- 
sional canonical system characterized by the inverse tem- 
perature /? = l/{kBT) and made up of identical particles 



of mass mo moving in the potential V{x). The stochas- 
tic element that appears in Eq. (^3), {B^, u > 0}, is a 
so-called standard Brownian bridge defined as follows: 
if {-Bti, M > 0} is a standard Brownian motion start- 
ing at zero, then the Brownian bridge is the stochas- 
tic process {Bu\Bi = 0, < m < 1} i.e., a Brow- 
nian motion conditioned on _Bi = S3|. In this pa- 
per, we shall reserve the symbol E to denote the ex- 
pected value (average value) of a certain random variable 
against the underlying probability measure of the Brow- 
nian bridge S". To complete the description of Eq 
we set Xr{u) — x + [x' ~ x)u (called the reference path), 
a = (/i^/3/too)^/^, and let pfp{x,x'] f3) denote the density 
matrix for a similar free particle. 

Rather than directly employing Eq. (Q, chemical 
physicist's arguments are usually constructed around the 
Trotter composition rule |0| that exploits the fact that 
{e~^^;/3 > 0} is a semigroup of operators on i^(IR), so 
that 



■P2H 



or, in coordinate representation. 



/^|g-(/3i+/32)ff|^'X ^ 



z){z\e 



-P-"\x') 



(2) 



(3) 
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By writing (3 — X)fc=i Pk, repeatedly applying the Trotter 
rule and choosing an adequate "short-time" approxima- 
tion, one ends up with a sequence of integrals on spaces of 
increasing dimension, converging to the density matrix as 
maxi<fc<„ j3k — > 0. Of course, this is much in the spirit of 
the original Feynman path integral approach. The meth- 
ods deduced by this technique are usually called Discrete 
Path Integral (DPI) methods (see Ref. [13 and the cited 
bibliography) . 
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It has become apparent that the Kag interpretation 
of the Feynman's formula may, in fact, offer a valuable 
starting point for the general construction finite dimen- 
sional approximations to the density matrix. This is so 
because the Brownian motion is a mathematically well 
understood object, for which various constructions are 
known. For example, the use of the Ito-Nisio theorem 
[Ts] has lead the present authors to the development of 
the random series path integral methods in a surprisingly 
general fashion [13 . This generality of the theory allows 
one to identify optimality criteria and eventually to an- 
swer questions as what the best representation is or how 
to modify the approach in order to improve the conver- 
gence. 

In this paper, we shall look at the relation between 
the Kag interpretation of the Feynman formula and dis- 
crete path integral methods. In the first part, we consider 
what we call the standard DPI methods. We shall again 
show the strength of the Kag approach, at least in terms 
of generality and mathematical interpretation of the for- 
mulae. In the second part, we explore the connection be- 
tween the random series technique, as particularized for 
the Levy-Ciesielski series representation of the Brownian 
bridge, and certain DPI results from the chemical liter- 
ature. While not the primary goal, we do obtain in an 
effortless manner a better way of implementing the latter 
results that features a "built in" fast sine-Fourier trans- 
form. Wc shall also derive the two basic modifications 
of the Levy-Ciesiclski path integral method (LCPI), the 
par tial averaging 13, 14J and the reweighted techniques 
[l3| . and establish their asymptotic law of convergence. 
We suggest that these results again emphasize the power 
of the Kag interpretation of the Feynman formula. By 
providing a central framework for discussion and analy- 
sis, the Kag approach significantly aids in characterizing 
the various methods and in establishing their intercon- 
nections, links that otherwise would be obscured by the 
multitude of possible representations. 



II. THE STANDARD DISCRETE PATH 
INTEGRAL METHOD 

A. Trotter-Suzuki approach. 

Our definition of the standard DPI method has to do 
with the particular short-time approximation that re- 
places the exact one in the Trotter product 



and 



rH 



rH 



(4) 



n-t-l terms 



We follow closely the arguments of Suzuki fis'l. The 
Hamiltonian of the system can be written as a sum be- 
tween the kinetic operator and the potential operator in 
the form H = K + V . The coordinate representations for 
the two operators are analytically known to be 

{x\e-P''\x')=pfj,{x,x'-p) 



{x\e-^^\x') ^e-^^''^^5{x' -x), 



respectively. It is therefore natural to consider short- 
time approximations that can be expressed by a finite 
composition of the above density matrices. The simplest 
example is the 2-term splitting formula 

^^f3{K+v) ^ e-^^e-'^^il + (5) 
which is of order 1. More generally, the order of a split- 



ting formula is said to be k if the relative error is 0{l3 
The motivation for this is that if 



fe+i 



then 



-P(K+V) 



1 + 



/5' 



fc+i 



(6) 



i.e., the error of the final n-tcrm Trotter product formula 
decays as fast as 1/nK The relation © was actually 
proved by Suzuki |l6j | in terms of operator norms for 
bounded operators A and B, but such an estimation also 
holds for K and V (which generally are unbounded op- 
erators; however, e~P^^~^^^ is bounded for most of the 
potentials of physical interest and for all positive (3). 
A better splitting is offered by the 3-term formula 



-P{K+V) 



.^0V-PK-^(3V 



6-3/3^1 + (7) 



or the one obtained by permuting V with K. These are 
of order 2 and go by the name of symmetrical trapezoidal 
Trotter short-time approximations [Tsl Il7j . More gener- 
ally, let us define a 21 + 1-term splitting formula of order 
k by the expression 

g-/3(if+y) ^ ^-aol3V^-bil3K^-ail3V 

___g-h,/3ifg-a,/3Vj^^0(^fe+l)j^ (8) 

Symmetry arguments suggest that for the optimum 21 + 
1-term splitting formula the sequences ao,...,a; and 
61,..., 5/ should be palindromic i.e., if the coefficients 
are read left to right, they form the same numerical se- 
quence as if they are read right to left. A look at the 
trapezoidal Trotter formula shows that this condition 
is natural, as one has little reason to believe that any- 
thing new can be achieved by considering some arbitrary 
^-ai3v^-l3K^-{i-a)i3v decomposition. In fact, with the 
help of the Campbell-Baker-Haussdorf-Dynkin formula 
[iSj, it can be shown that this more general expression 
is an order 2 splitting only if a = 1/2 and that it is an 
order 1 splitting otherwise. More generally, since the op- 
erator e~P^ is Hermitian, it is natural and, as argued by 
De Raedt and De Raedt |l7j, optimal to approximate it 
by a sequence of Hermitian operators. It is straightfor- 
ward to see that the n-order Trotter product Q is Her- 
mitian if and only if the short-time approximation ||SJ) is 



3 



Hermitian. In turn, this requires the palindromicity of 
the {oi} and {bi} sequences. 

It is not difficult to see that H J^i'^i — P ^^d — 
q, then the n-order Trotter formula Q converges to 
exp[— /3(gif + On the other hand, the equality 

holds for arbitrary potentials V{x) if and only if q = 1 
and p — 1. Therefore, the additional constraints = 
1 and J2i = 1 must be enforced upon the sequences 
{a,} and {b,}. 

We considered this more general problem with the 
hope that by using a more advanced splitting one may 
improve the asymptotic order of the Trotter product for- 
mula. Now, there is one more restriction that we have to 
place on the sequence ao, 6i, ai, . . ., namely it should be 
made up of real and positive numbers only. Otherwise, 
the short-time approximations are either ill-defined or, by 
Trotter composition, generate algorithms that are numer- 
ically unstable at low temperatures. Unfortunately, the 
following theorem of Suzuki (see Theorem 3 of Ref. ITsj) 
says that 

Theorem 1 (Suzuki nonexistance theorem) There 
are no finite length splitting formulae ^ of order 3 or 
more such that the coefficients ao,bi, ai, . . . are all real 
and positive. 

This negative result shows that more general splitting 
formulas do not produce short-time approximations ca- 
pable of improving upon the trapezoidal Trotter result, 
at least as far as the asymptotic order of the Trotter 
product rule is concerned. However, the product rule 
which uses equally spaced time slices, does not pro- 
vide the most general standard DPI expression. In the 
next section, we shall argue that this most general ex- 
pression is of the form given by Eq. (jSJ), for which the 
Suzuki nonexistance theorem does not apply. 

B. Direct quadrature of the Feynman-Kag formula. 

Let us notice that the form of the equation |(SJ) is in- 
variant under the Trotter composition Q and so it can 

I 



One may use the above joint distribution density to com- 
pute the expectations of the functionals of the Brownian 



be regarded as the most general standard DPI approxi- 
mation to the density matrix provided that we can give 
a recipe for choosing the sequences ao, 6i, ai, . . . , 6;, a; in 
such a way that the correct result is recovered in the limit 
I —^ oo. While in the Trotter-Suzuki approach this may 
seem a daunting task, the problem has an easy solution 
by means of the Kag interpretation of the Feynman for- 
mula. In this section, we shall derive a more general ex- 
pression for the standard Discrete Path Integral method 
simply by replacing the monodimensional integral over 
u in Eq. with an approximate quadrature sum and 
then using the definition of the Brownian bridge to com- 
pute the expectation of the resulting functional. Given 
the Suzuki nonexistence theorem, it is hard to believe 
that one may eventually devise a standard DPI method 
with asymptotic convergence 0{l/n^) or better. How- 
ever, before one starts to investigate the validity of this 
conjecture, one needs a more general statement of the 
standard DPI method. 

For obvious reasons, the random process W^^,{u) = 
Xr{u) + (jB^ is called a Brownian bridge of variance 
and end points {x,x'). The Feynman-Kag formula can 
be expressed in terms of the new process in the form 

^^f^^ = Eexp (-/? rv\wiAu)]du] . (9) 
Pfp{x,x';(i) I Jo ^ -I J 

A quite important property of the Brownian bridge 
W^^,{u) is the joint distribution of the variables 
■ • ■ ' ^:^',x'('^n) fo'" ^ given partitioning < 
ui < . . . < u„ < 1 of the interval [0, 1]. Let us set 



Pt{x) = 




and notice that pfp{x, x' ; f3) = p^2 (x' — x). From the very 
definition of the Brownian motion [T9l |. the aforemen- 
tioned joint distribution can be straightforwardly shown 
to be 



(10) 



I 

bridge which are of the form 



P {Wx,x'i'>J-l) ^ [xi,Xi + dxi], . . . , W^.j.,{Un) e [Xn,Xn + dXn]} 
^Pa^uA^l - x)Pa^(u2-ui){x2 - Xi) . . . ^^^2 („„_„„_ ^ ) (a;„ - Xn~l) 
XP<t2(i_„^)(x' - Xn)/Pa^ix' - x) dxi . . .dXn- 
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p,2{x' -x)E{f[W:,,{u,),...,W:^,{u^)]}= dx,... dXnf{xi,...,Xn) 

xpcr2„i(a;i - a;)pcr2(„2-„i)(a;2 - xi) . . .po-2(u„-u„_i)(2:n - Xn~i)Pcr^{i-.u„)ix' - a;„), 



(11) 



where f{xi,...,Xn) is some integrable n-dimensional 
function. As a direct application of Eq. Hll|l . consider 
a quadrature scheme on the interval [0, 1] specified by 
the points = uq < ui < . . . < Un < Un+i = 1 and 
the corresponding nonnegative weights wq,wi, . . . , Wn+i- 
Replacing the monodimensional integral in the Feynman- 
Kag formula by its quadrature form, we obtain an 
approximation to the density matrix of the form 



o^^\x,x'-l3) 



n+l 



i=0 



du 



(12) 

The expectation value of this formula can be exactly re- 
duced to a finite dimensional integral with the help of 
the formula itTT)) . We call Eq. ifT^ the standard Discrete 
Path Integral (DPI) method and we expect it to converge 
to the correct result for all continuous and bounded from 
below potentials V{x). In this respect, remember that 
with probability one the Brownian paths are continu- 
ous and therefore so is as a function of u. 
Also remember that by definition a quadrature scheme 
on [0, 1] is constructed so that it eventually integrates all 
continuous functions on [0, 1]. 

Formula H12|) can indeed be formally deduced start- 
ing with the Trotter composition rule (|2Jl and a carefully 
chosen sequence of short-time approximations. More pre- 
cisely, for z = 0, 1 . . . , n, define 9i = Mi+i — Ui. Then the 
equation l(T^ is nothing else but the Trotter product 



^-wo0V ^-BoPK ^-wiPV ^-eil3K 



(13) 



which is of course of the type given by the formula (jS)). 
Finally, let us notice that we always have Oi — 1. We 
also require that ~ 1- 1^ fact, one is not interested 

in relaxing these equalities because they are the necessary 
and sufficient conditions to obtain the exact free particle 
density matrix at all levels of approximation. The reader 
can directly verify this fact by assuming that the poten- 
tial V(x) in Eq. is constant but not zero. Moreover, 
the additional restriction of the integration schemes to 
those for which the sequences Oi and Wi are palindromic 
is justified by the requirement that the approximate den- 
sity matrices be Hermitian. 

To summarize, the advantage of the equation H12I) is 
the novel interpretation for the sequences 6i — Ui+i — Ui 
and Wi, leading us to a more general convergence prob- 
lem: what is the best convergence order for the stan- 
dard DPI approach H12|l and for what types of quadrature 
schemes is it attained? As we suggested in the beginning 



of this subsection, it is very plausible that the answer 
to the above question is two and is attained for almost 
all "sensible" quadrature schemes. We illustrate this by 
studying the convergence of the diagonal matrix element 
p[0-(j) = {Q\e-I^"\Q) of an harmonic oscillator for the fol- 
lowing quadrature techniques: the trapezoidal rule (TT) 
and the Gauss-Legendre method (GL) '2^1 . In both cases 
the condition Wi — 1 and the palindromicity of the 
sequences 9i and Wi are respected. We leave it for the 
reader to show that if the trapezoidal rule is used for 
integration, then one recovers the classical trapezoidal 
Trotter formula. 

If Mt stands for any of the methods studied and if 
aut represents the convergence order of the correspond- 
ing matrix element p*^*(0; /3), then the convergence con- 
stant is defined by 



CMt = lim n"«*[p(0;/3) 



Mt 
rn 



(0;/?)]. 



The above relation can be cast in the more intuitive but 
equivalent form 



p(0;/3)«pr(0;/?) 



CMt 



with an appropriate definition of the symbol w. These 
convergence orders and convergence constants can be 
evaluated numerically as follows. For each method, we 
compute 



{r? - 1/4) log 



P4T2(0;/3)-pe+2(0;/?) 

P4'^V2(0;/3)-p(0;/3) 



where P4^*+2(0;/3) represents the DPI approximation of 
order An + 2 for the method Mt. The evaluation of the 
matrix elements /54;^+2(0; /3) is discussed in Appendix B. 
Then, as argued in Ref. ifl^ . a^^* as a function of n is 
asymptotically a straight line, whose slope gives the con- 
vergence order. Therefore, aMt — lim„^oo ct^+i ~ ct-n* ■ 
As to the convergence constants, they can be evaluated 
by studying the asymptotic slopes of 



Mt 



{An + 2)"-'(n + 1/2) K*,2(0;/3) - p(0;/3)] , 



once aMt is known. The computations were performed 
in atomic units for a particle of mass mg = 1 and for 
the harmonic oscillator V{x) = x^/2. The inverse tem- 
perature was (3 — 10. As shown in Fig. ^ the asymp- 
totic convergence order of both methods is 2. The con- 
vergence constants are found to be ctt — 0.103 and 
cgl = 0.127, respectively. One notices that at the 
temperature (3 = 10, the trapezoidal Trotter method 
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FIG. 1: The current slopes q*+i — a*^* for the trapezoidal 
rule (TT) and for the Gauss-Legendre quadrature method 
(GL) are shown here to converge to the same value of 2. 



is slightly faster. However for /3 = 1, one computes 
ctt = 0.033 and cql = 0.005, which indicates that for 
this temperature the Gauss Legendre method is faster. 

The conclusion we draw from this analysis is that "bet- 
ter" integration schemes do not necessarily improve upon 
the convergence of the standard DPI methods. Why is 
this so? Again the Kag interpretation of the Feynman 
formula gives us an explanation which is not obvious from 
the Trotter composition rule. A famous theorem due to 
Paley, Wiener, and Zygmund [23| says that with prob- 
ability one the paths of the Brownian motion are con- 
tinuous but not differentiable at every point. Therefore, 
V[xr{u) -f ^ ^ function of u is not differentiable even 
if the potential V is. As emphasized by Press et al. [20l |. 
higher order quadrature schemes do not automatically 
translate into better convergence, unless the integrand is 
well behaved. In our case, there is a limit upon the rate 
of convergence of the quadrature schemes which is set by 



the properties of the Brownian motion paths rather than 
by the properties of the potential, provided that the lat- 
ter has a continuous first order derivative. In conclusion, 
one expects that there is an intrinsic limit for the conver- 
gence order of the standard DPI methods. Moreover, the 
Suzuki nonexistence theorem predicts that none of the 
classical quadrature formulas for equally spaced abscis- 
sas (e.g. Simpson's rules, Bode's rule, etc.) are going to 
improve upon the asymptotical convergence of the trape- 
zoidal rule. This is strong evidence that the intrinsic limit 
for the convergence order of the standard DPI methods 
is 2. 



C. Kinetic energy diagonalization for the standard 
DPI technique 

In practical applications, it is generally difficult to work 
directly with the formula because this involves a 
correlated Gaussian multidimensional distribution. As 
shown by Butler and Friedman 22], this correlated dis- 
tribution can be replaced with an independent one by 
simple algebraic manipulations. Later, Coalson [2^ used 
a similar technique in order to demonstrate, on an in- 
tuitive basis, the relation between the discrete and the 
Fourier path integral methods. The two approaches men- 
tioned above are technically different and in fact there are 
an infinite number of such transformations. As we shall 
see in this section, they are related by simple orthogonal 
transformations and in Section III we shall propose a new 
approach, which allows for faster numerical implementa- 
tions. 

We begin by performing a coordinate transformation 
so as to diagonalize the positive definite quadratic form 
associated with the kinetic operator. More precisely, 
let us introduce the transformation of coordinates z„ — 
{un)]/(J. By using the condition X)*^* = 1; is 
straightforward to show that the formula 1)11(1 becomes 



dzi . . . y dz„ /[Zicr -I- Xr{ui), . . . , Z„CT -I- Xr{Un)] 
Xpg„{zi)pgjz2 ~ Zi).. .pg^_^{Zn - Z„_l)p6l„ (^^n), 
I 



or, in an even more compact notation, 

E{f[WlAni),---,Wl,>{un)]} 

1 



dzi 



v/(27r)"det(A)-i 

x/[ziO- -I- Xriui), . . . , Z„a + XriUn)], 



exp 



^-z'^Ai 



where the matrix A is the n-dimensional tridiagonal ma- 
trix defined by Ai i — l/9i + l/6'i_i for 1 < i < n and 
A,^i+i = Ai+i^i = -l/9i for 1 < i < n - 1. 



By construction, the matrix A is symmetric 
and positive definite [otherwise, the integrability of 
exp {~z^ Az/2) would be violated] and can be diagonal- 
ized by an orthogonal matrix S. Defining the new coor- 
dinates y = S'^z, and letting {A;; 1 < i < n} be the set 
of the n real and (strictly) positive eigenvalues of A^ we 
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have 

nfmA^i),...,w:,,,iu„)]} = 

X / dyi 



-I 1/2 



/ dyi . . . / dy„ exp ( V Aiy,2 ) 



i=l i=i 



1/2 

Finally, setting Ui — yi, one ends up with 
E{f[Wl,,{u,), W^,,,{u,,)]} = / dai . 



Xr{ui) 



da„ 
(14) 



+^ J TT72 , • ■ • , a;^ (w„) + cr ^ S'n.j 



,1/2' ■ ■ ■ '•"'^^"'"^ , ,/2 

This formula is advantageous for numerical applications 
because the integration is performed over independent 
identically distributed Gaussian distributions. Given a 
quadrature scheme, one diagonalizes the tridiagonal ma- 
trix A and tabulates the values of Sij and A^. For the 
case of equally spaced time slices, the eigenvectors and 
eigenvalues of the matrix A are known analytically: 



2 . / IJTT 

■ sin 



n+1 



and 



Xi = 4(n + l)sin^ 



n+1 



_2(n+ 1) 



1 < i, j < n 



1 < i < n, 



respectively. 

Similar to the invariance of the Brownian bridge at a 
change of basis as shown by the Ito-Nisio theorem, the 
formula p4|l is invariant to arbitrary orthogonal transfor- 
mations of the vectors a = (ai, . . . , a„). Indeed, let Q be 
an arbitrary n-dimensional orthogonal matrix, and con- 
sider the coordinate transformation a' = Q^a. Notice 
that J2i = 12i ^'1 s-i^d define the matrix 



-''J ~ ^l 



k=l \- 



/2 



Then a little algebra shows that 

x(2^)-"/^xp(^-iXJa,^^/[x.(ui) 

n n 

i=i i=i 



(15) 

da„ 
(16) 



Because of the additional degrees of freedom, the last 
formula is more useful in practical applications than the 
transformation A good part of the computational 

time is spent with the evaluation of the current paths. 
For a monodimensional system, one usually needs a num- 
ber of operations proportional to in order to com- 
pute the vector T a by matrix multiplication. However, if 
equally spaced time slices are used, the n elements of the 
form J2]Li -Sijaj/A^^ from Eq. ^ can be computed 
by fast sine- Fourier transform in a number of operations 
proportional to nlog2{n), provided that n = 2*^ — 1 with 
A: > 1 [Hi^. Equivalently, one may say that there 
must be some orthogonal matrix Q such that the asso- 
ciated matrix T defined by the relation H15|l is a sparse 
matrix with at most k nonvanishing elements on any line. 
Therefore, the evaluation of the elements Ta by direct 
matrix multiplication requires only 0{k ■ n) operations. 
In this paper, we shall directly find such a matrix T by 
means of the Levy-Ciesielski representation of the Brow- 
nian bridge, which is discussed in the next section [see 
formula (|77|) ]. 



III. THE LEVY-CIESIELSKI 
REPRESENTATION OF THE FEYNMAN-KAg 
FORMULA 

As we discussed in the previous section, the transfor- 
mation (|14|l was utilized by Coalson in order to establish 
a connection between the discrete path integral meth- 
ods and the Wiener- Fourier path integral technique 2^^ . 
However, strictly speaking the Wiener-Fourier sequence 
of approximations is not equivalent to any discretization 
scheme. That is, for any n, there is no sequence of short- 
time approximations which by Trotter composition would 
generate the nth order Wiener-Fourier approximation. A 
more precise statement of this assertion is given at the 
end of Section III.B. Then, a natural question arises: Is 
there any random series for which at least a particular 
subsequence can be thought of as a DPI method? The 
answer is positive and is furnished by the Levy-Ciesielski 
random series construction of the Brownian bridge. 

In this section, we shall specialize the general the- 
ory of the random series representation of the Feynman- 
KaQ formula ^3 for the particular case of the Levy- 
Ciesielski representation of the Brownian motion. The 
respective method will be designated by the acronym 
LCPI. We shall also derive the three associated meth- 
ods: the primitive LCPI, the partial averaging LCPI, 
and the reweighted LCPI. Moreover, with the help of the 
Levy-Ciesielski series representation, we shall prove the 
Trotter product rule for the case n — 2'' — 1 and for this 
subsequence, we shall show that each of the above mod- 
ifications of the LCPI method can be interpreted as the 
n-order Trotter product of some appropriate short-time 
approximations. In doing so, we establish a direct con- 
nection between the discrete and the random series path 
integral techniques. As a practical application, we shall 
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obtain a sparse matrix T of the form (|15|l which requires 
only 0{k ■ n) operations to compute the vector Ta by 
matrix multiphcation. 

A. The Levy-Ciesielski path integral method. 

Some of the arguments we use in the following in- 
troduction to the Levy-Ciesielski re pre sentation of the 
Brownian bridge can be found in Ref-^S*. For k — 1,2, . . . 



0.50H 




u 



FIG. 2: A plot of the renormalized Schauder functions for 
the layers k — 1,2, and 3 showing the pyramidal structure. 

are called the Schauder functions. As McKean puts it 
|2^, the Schauder functions are "little tents," which can 
be obtained one from the other by dilatations and trans- 
lations. In modern terminology, this has to do with the 
fact that the original Haar wavelet basis is a multiresolu- 
tion analysis of i^([0, 1]) organized in "layers" indexed by 
k j2^. If we disregard the factor 2^*^"^)/^, the Schauder 
functions make up a pyramidal structure as shown in 
Fig.H 

Let {ofcjsfc l,2,...;j ^ 1, 2, . . . , 2*^-1} be i.i.d. 
standard normal variables and define YQ{u,a) = and 

Yk{u,a) = ^ ak,jFk^j{u). 
Then by Ito-Nisio theorem, 

oo 

B°{a) =Y,Yk{u,-a) (19) 

fc=i 



and j = 1, 2, . . . , 2*^ ^, the Haar function fk j is defined 
by 

r 2('=-iV2^ te [(^-l)/2^^/2'=] 

fkAt) = \ -2('=-i)/2, t G [V2^ {I + l)/2'=] (17) 
[ 0, elsewhere, 

where I = 2j — 1. Together with /o = 1, these func- 
tions make up a complete orthonormal basis in L^{[0, 1]). 
Their primitives 



(18) 

I 

is equal in distribution to a standard Brownian bridge 
and the convergence of the right hand side random series 
is uniform almost surely. 

Let us now define the primitive, partial averaged, and 
reweighted LCPI methods, which are the standard tech- 
niques that can be derived from a series representation 
[13. They will be denoted in the following discussion 
by the acronyms Pr-LCPI, PA-LCPI, and RW-LCPI, re- 
spectively. The nth order Pr-LCPI term is obtained by 
approximating the Brownian bridge by the n-dimensional 
process 

k j 

Suia) = ^ri(M,a) -H^afc+i,iFfc+i,i(w), 
(=1 (=1 

where k and j are the unique numbers such that n — 
2*^ + j - 1, with A: > and 1 < j < 2*^. However, it 
appears natural to utilize only the subsequence of the 
form n = 2*^ — 1 with fc > 0, corresponding to k complete 
layers and from now on we shall restrict our attention to 
this subsequence, for which 

k 

S:{a) ^Y.^iiu,a). (20) 

/=i 

Using the notation introduced in Ref. 0, we denote 
the tail of the series ((T^ by 

00 

l=k+l 

To define the PA-LCPI method, besides the sum (|20|l . we 
need to evaluate 

00 2 
i=k+i j=i 



2(k-i)/2[t_ ^l_iy2% te [(^-l)/2^V2''] 

Fk^jit) = { 2(^-i)/2[(; + i)/2'=-t], t e [^/2^ (? + l)/2'=] 

0, elsewhere 
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This quantity must be computed explicitly because it 
enters the final PA-LCPI formula by means of the "effec- 
tive" potential 



Vu.nix) 



: exp 



V{x + z)dz. 



For more information, the reader is referred to the Sec- 
tion III of Ref. m For all I > k + 1 and 1 < j < 2\ 
the functions Fij are zero on the points Up — with 
p = 0, 1, ... 2'^. Let us define the support of the function 
F;j as the set supp(F;j) = {u e [0,1] : Fij{u) ^ 0}. 
Moreover, for 1 < p < 2^ let Ip = {{IJ) : l>k + l,l< 
j < 2', supp(F;.j) C [up^i,Up]} and define 

Wp{u,a)^ ^ aijFi^jiu). 

Then a little thought and the use of the Ito-Nisio theorem 
show that Wp{u, a) is a Brownian bridge on the interval 
[up_i,Up] of variance 1/2'^. In addition, if pi P2, then 
the Brownian bridges Wp-^ {u, a) and Wp^ {u, a) are inde- 
pendent because they are functions of the independent 
Gaussian random variables {ai,j} with G Ip^ and 

G /p2, respectively and the sets of indexes Ip-^ and 
/p2 are disjoint. It is convenient to denote by Ep the con- 
ditional expectation over the random variables aij with 

G Ip. Then, we have 



Ref.m the main idea is to simulate the effect of the par- 
tial averaging method by replacing the tail series i?"(a) 
in the full series expansion by a collection of random 
variables {i?"(&i, • • • , 6„+g)}o<u<i defined over an n + q 
dimensional probability space (g is a small integer which 
does not depend upon n). We ask that (i) the variance 
at the point u of • • ,6„+g), denoted by r^(u), 

be as close as possible to r^(M) and (ii) the variables 
5"(ai, • • • , a„) and • • • , 6n+g) be independent and 

their sum have a joint distribution as close to a Brown- 
ian bridge as possible. One candidate for our approach 
is • • • , 6„) = X]p=i hpVlpiu), where 6i, • • • , 5„ are 

i.i.d. standard normal random variables. Condition (ii) 
above is realized in the Ito-Nisio theorem by insuring 
that the collection {i^/.j (u), with 1 < ? < fc, 

1 < j < 2'^^, and \<p<n + q\s orthogonal and 
we shall look for such a collection. Here, LL!p{u) is the 
derivative of r2p(u) and it is not required to be normal- 
ized. 

As opposed to the Wiener-Fourier series, for the Levy- 
Ciesielski series it is possible to enforce the condition 
(i) exactly. The analysis done for the partial averaging 
method showed that we can represent r^(M) as the sum 
of 2*^ = n + 1 functions of disjoint support and which 
can be obtained one from the other by translation. In- 
tuitively, we must set g = 1 and replace the 2^ Haar 
functions making up the fc + 1 layer by 



and 

nBZ(af]^Y.^[Wp{uraf]^Y.^p[Wp(u,a) 
However, one computes 



(21) 



2 / \ m ni^ / -\2^ \ uil - 2''u), < W < 2 



0, 



otherwise 



(22) 



and then by translation 

^lp{u) ^ ¥.p[Wp{u,af] = ^l^[u~{p-l)/2\ (23) 

Clearly, the functions 7^p(u) have disjoint support. Fi- 
nally, Eq. becomes 



r2(«) = a2^7^.p(«), 
p=i 



(24) 



which concludes the definition of the PA-LCPI method. 

The reweighted technique is yet another way of im- 
proving the convergence of the primitive method. It has 
the advantage that it does not require the evaluation of 
the Gaussian transform of the potential. As discussed in 



Wp(u) = ^7n[u-(p-l)/2'']. 

It is easy to notice that the functions ujp{u) are orthog- 
onal among themselves because they have disjoint sup- 
port. Moreover, it is not difficult to see that the Haar 
functions fi,j{u) are constant on the intervals [up_i,Up] 
for all I < k and therefore, they are orthogonal on the 
ujp{u) functions because 



Wp(u)du = / Wp(u)du = 7„,i(l/2 ) - 7„4(0) = 0. 

Jxp-i 



In consequence, the n-order RW-LCPI approximation 
uses the series 



k 

E 

1=1 



Yi{u,a) 



2" 



(25) 



for its implementation. [It is customary to define the 
approximation order by the dimensionality of the under- 
lying probability space. We shall not apply this rule in 
the present paper in order to keep the unity of the ex- 
position. The squares of the functions 7n,p(ii) are given 
by the relation (|23|l .] A look at Fig O shows that the n- 
order RW-LCPI method is identical to the 2n + 1 order 
Pr-LCPI, except for the replacement of the last layer of 
fimctions Fk+i.p{u) with "fn,p{u). 
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0.50 



0.25 



0.00^ 




where the short-time approximations p^* (x„, a;'; /?) are 
defined as follows: 



exp 



-P J V[x + {x' - x)u]du } , 



pE\x,x'-f3) r /-i 

nK^^'^^\~"l V u.q[x + [x - x)u\au 

Pfp{x,x';f3) { Jo 



FIG. 3: A plot of the functions used in the reweighted LCPI 
technique of order 2. Note that the "httle tents" of the layer 
k = 3 were replaced by "little domes." 



B. Properties of the Levy-Ciesielski path integral 
method. 



As announced in the beginning of the section, the LCPI 
method for n = 2*^ — 1 is virtually a reformulation of the 
Discrete Path Integral method with appropriate short- 
time approximations. This is shown by the following re- 
sult. 

Theorem 2 If n — 2*^ — 1, then the following relations 
hold true: 

1. The Trotter theorem 

p{x, x'\ /?) = J dxi.. . J dXn p Xi] Jl^^ 

( /. p 

...pi X'fi , X , 



n+1 



2. If Mt stands for any of the LCPI methods, then 



p^'\x,x';P)^ dxi... da;„p*" 



n+1 



Mt I I. P 



n -I- 1 



prix,x';(3) 
Pfp{x,x';P) 



dz(27r) 



-1/2 -.V2 



xexp'l— /3 J V[x + {x ~ x)u + zTo{u)]du 
respectively. 



Proof. We only prove the first point of the theorem. 
It is not difficult to see that the Trotter theorem is 
in fact a part of the latter case with p^'^{x,x'](3) = 
p^'*{x,x'; P) = p{x,x';P). As such, the second point 
follows by arguments similar to the first one and is left 
to the reader. 

Let us remember that for alH > fc -I- 1 and 1 < j < 2' , 
the functions Fi,j{u) are zero on the points Up — p/2'^ 
with p = 0, . . .2*^. This means that the joint distribution 
of the Brownian bridge at these points is uniquely deter- 
mined by the series (|20|) . For this proof, it is important 
to notice that the inverse result is also true: knowledge 
of the joint distribution of the points Up with p = 0, . . . 2*^ 
uniquely determines the series H20() because the latter is 
linear on the intervals [wp-i, Up]. It follows that the vari- 
ables are independent of the displaced and rescaled 
Brownian bridges Wp{u,a). Using this information to- 
gether with the joint distribution density for the random 
variables Xr{up) -\- crB'^^, which is given by the formula 
as shown in the previous section, one computes 



p<,2(a;' -a;)Eexp<{ -/3 / V Xr{u) + aB° dw J> = p„2{x' - x)Ecxp I -p"^ f V 



Xr{u) + aB^ 



du 



dxi ■ I dxn W P_a±_ {xp+i - Xp)¥.p exp I -/3 



p=0 



V 



{u ~ Up) 
(Up+l - Up) 



{xp+i - Xp) + aWp{u, a) 



du 



dxi 



p=0 



1 • • • / dxn TT ( p^ {xp+1 - Xp)E exp 



P 



n+l Jo 



V 



€p + u{xp+i - Xp) + B,^(a) 



du 
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which proves the first point of the theorem. The latter 
follows by a similar line of thought: for instance, the 
result for the primitive method is obtained by setting 
Wp{u,a) = in the previous formula. □ 

From the theoretical point of view, the importance of 
the Theorem 2 consists of the fact that it establishes 
a direct connection between the random series and the 
discrete path integral techniques, even if only for the 
n = 2*^ — 1 subsequence. As such, we notice that the prim- 
itive result was employed by Makri and Miller |27ll2q and 
by Mielke and Truhlar [0| as the ZOP-DPI method. The 
latter authors found that the asymptotic convergence of 
the method was 0(l/n). This result is in good agreement 
with the present analysis because the primitive Levy- 
Ciesielski method cannot exceed the convergence rate of 
the most rapidly convergent series, the Wiener-Fourier 
series, which behaves asymptotically as 0(l/n) [l3j |. 

The partial averaging result is not new either. The DPI 
formulation was used by Kole and De Raedt to treat 
systems with negative coulombic singularities, for which 
the non-averaged methods are known to be ill-behaved. 
However, Kole and Dc Raedt were not aware of the fact 
that they were using the partial averaging method in a 
special setting and regarded their approach as a separate 
method. It has been shown in a mathematically rigorous 
way that the partial averaging method is convergent for 
such potentials at least as far as the pointwise density 
matrix, the partition function, and related integral ex- 
pressions are concerned 30] , for all series representations 
of the Feynman-Kag formula. Therefore, the Theorem 2 
can be used to give a mathematically rigorous proof of 
the Kole and De Raedt result, which conversely can be 
thought of as an argument demonstrating the desirable 
properties of the partial averaging strategy. 

In a related situation, the partial averaging method as 
specialized for the Wiener-Fourier series representation 
was used to treat the polaron problem by Alexandrou, 



Fleischer, and Rosenfelder [s^l- Later, the DPI formu- 
lation of the partial averaging technique was applied by 
Titantah, Pierleoni, and Ciuchi l^^l for the same polaron 
problem and regarded once again as a separate tech- 
nique. We hope we have convinced the reader that given 
the multitude of series representations that may enter 
the Feynman-Kag formula, there are an infinite number 
of ways in which the partial averaging idea can be im- 
plemented. The Wiener-Fourier and the Levy-Ciesielski 
series representations as well as the related DPI imple- 
mentation are only some instances (although perhaps the 
most important ones). 

In Section II. C, we promised that we would find a quick 
way to compute the current paths for the standard DPI 
methods by means of the Levy-Ciesielski series represen- 
tation. For the LCPI formulation, it is straightforward 
to notice that the computational time necessary to com- 
pute the current path at a point u is proportional to 
k = log2(n -f 1) for the Pr-LCPI and PA-LCPI methods 
and H-log2(n-f 1) for the RW-LCPI method, respectively. 
This is so because given a point u, the only Schauder 
function from the layer / that is non-zero at the point 
u is Fij{u) with j = [2^~^u] + 1, where [x] denotes the 
integral part of x. For the RW-LCPI method, we have in 
addition that the only function 7„ (u) which is non-zero 
at the point u is the one with j — [2*^^] -I- 1. In fact, 
going back to the proof of Theorem 2, we remember that 
the joint distribution of the points Up with p — 0, ... 2*^ 
uniquely determines the series H20|l because, in a more 
mathematical notation, we have 

k 

Sl{a)^Bl^{a)^J2^i{up,d), V 1 < p < n. (26) 
1=1 

Equation 126|1 allows us to write the following special 
form for Eq. (fT?)|l : 



/ dai . . . / da„ (27r) "^^ exp ( XI 



X/ 



1=1 i=l 
fc 



-P';,[2'-iMi] + l("l)ai,[2'-i«i]-|-l, ■ • ■ ,Xr{Un) + O" ^ ^i,[2'-i«„] + l ("«)«;, [2'-iti„] + l 

I 



1=1 



1=1 



This proves that the standard DPI method can be im- 
plemented so that the number of operations necessary to 
compute the current paths is 0{k ■ n). 

As we said at the beginning of the Section III, as op- 
posed to the n = 2'' — ! subsequence of the Levy-Ciesielski 
representation, no subsequence of the Wiener-Fourier 
representation can be rationalized as a DPI method. The 



precise meaning of this is that if 



2E 

fc=i 



flfc- 



sin(A:7ru) 



(27) 
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then there is no sequence = 7io < ui, . . . . u„ < 
such that 



Convergence of the PA-DPI and of the RW-DPI 
methods 



E{.f[Wl,,{m),...,W^^,,{ur.)]} 
= E{f[xr{ui)+aS'^^{a),...,Xr{un)+aS:ja)]} (28) 

for all functions /(xi, . . . , Indeed, remembering 

W^^,{u) = Xr{u) + (jB^ and choosing f{x) — [x — 
Xr{up)]/ a"^ for some interior point < < 1, one com- 
putes 

nmiA^i). . . . , wiAun)]} = E«)^ = - uj,). 

On the other hand, 

nn^Au,) + oSZ^ (a), . . . , Xr{un) + a S^,^ (a)]} 



Clearly, the equality 



fc2 

cannot hold because 



Up{l - Up) 



■E 

fc=i 



SIV? [klTUp) 



■ E 



sin^(A:7rMp) 



does not vanish on the interval (0, 1). To prove this, it is 
enough to notice that the zeros of sin^[(n + l)7rup] and 
sin^[(n + 2)-KUp\ are strictly interlaced. 

The following theorem, whose proof is left to the 
reader, provides the necessary and sufficient conditions 
for an n-order term of an arbitrary series to admit a par- 
ticular m-order DPI representation. 



Theorem 3 Let 



u{l — u) 



k=l 



and let — uq < ui < . . . < u,„ < u„i^i = 1. Then 
= E{f[xr{ui) + aSZ, (a), • . ■ , -f aSl^ (a)]} 



for all f : 
1 < p < m. 



if and only if r^(up) = for all 



The fact that the Wiener-Fourier representation can- 
not be rationalized as a DPI method should not be sur- 
prising. Indeed, we presented enough evidence in Ref.'lS' 
to support the idea that the convergence of the partial av- 
eraging and the reweighted Wiener-Fourier path integral 
methods is 0{l/n^) for sufficiently smooth potentials. 
On the other hand, the analysis performed in Section II 
suggests that we cannot expect an asymptotic conver- 
gence of the DPI methods better than 0{l/n^). In fact, 
as we will show in the next subsection, the n = 2^^ — 1 
subsequence of the PA-LCPI and RW-LCPI methods can 
have at most 0{l/n'^) asymptotic convergence. 



In this subsection, we shall study the convergence of 
the Trotter product formulae having as short-time ap- 
proximations the partial averaging and the reweighted 
zero order formulae given in Theorem 2. It is natural to 
call this methods the PA-DPI and the RW-DPI methods, 
respectively. In particular, by virtue of Theorem 2, we 
obtain the asymptotic rates of convergence for the subse- 
quences n = 2'"' — 1 of the corresponding LCPI methods. 
To anticipate, the convergence of the partition function 
and of the density matrix will be shown to be 0{l/n^) for 
both methods if the potential is smooth enough. More 
precisely, we limit our discussion to the class of potentials 
introduced in Ref. lsoL which are the Kato-class potentials 
[sij having finite Gaussian transform. In this section, a 
potential is called smooth if it lies in the local Sobolev 
space W^];^{R'^) and if the squares of the potentials and 
of the first order derivatives have finite Gaussian trans- 
form. We remind the reader that the local Sobolev space 
Wl;^{R'^) is made up of all Lf^^'^) functions whose first 
order distributional derivatives are also Lf^^{U.'^) func- 
tions i.e., 



D 



y(.T)2 +^ |ay(a;)/5a:. 



da; < 



for all bounded domains D C M.'^. We warn the reader 
that the 0{l/n'^) convergence of the density matrix and 
of the partition function for this class of potentials does 
not automatically imply similar convergence for the en- 
ergy estimators, for which additional restrictions upon 
the class of potentials might be necessary. 

To simplify the notation, we prove the convergence 
results for the monodimensional case and only state 
the multidimensional analogues. Let us start with the 
asymptotic convergence of the partial averaging method. 
If we set 



Uix,x',P;a) = / V + crB°(a) du, (29) 

Jo L 



a little algebra shows that 



p{x, x';f3)- p^^{x, x'; /?) = p{x, x'; (3) - p^^{x, x' ■ (3) 



Po 



(x,x';/3)E|f 



'l3[U{x.,x' ,P;a)-EU{x.,x' .p-.a)] 



The first equality follows from the fact that zero order 
PA density matrix is always smaller than the true density 
matrix, according to equation (18) of Ref. However, 
for P small, we can expand the exponential in a Taylor 
series in order to establish the order of the short-time 
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approximation. We have: 



In particular, the inequahties 



p{x, x';f3)- p^^(.T, x';(3)^ p^^{x, x'; f3) 

X I [U{x, x', (3;a)-E U{x, x', /?; a)f (30) 

+ [U{x, x', f3;a)-E U(x, x' , (3; a)]'^ + 0{I3^] 

Notice that the term of order one in the Taylor expansion 
cancels, so the asymptotic behavior is dictated by the 
variance of the function U(x,x' , fi;a). However, looking 
at the expression (|29|) . we see that this variance must 
also decay to zero as /? ^ 0, because cr ^ 0. The same 
is true for the third order moment and a gradient expan- 
sion similar to the one performed in Appendix A for the 
variance of the function U{x,x' , (3;d) shows that 



p^''{x,x';(3)<p{x,x';(3) 



<p^^{x,x';(3) 



l + ^T2ix,x';p)+OiP^) 



show that the zero order partial averaging formula is of 
convergence order 2. Therefore, the assertion of Makri 
and Miller [23| that Pq^{x, x' \ 0) is not an order 2 short- 
time approximation is wrong. Also, notice that 



||-E [U{x,x',P;a) -EU{x,x',P;d)] 

decays to zero as fast as 0{f3'^-^). 

We shall be more careful in establishing a proper bound 
on the variance of the function U{x, x' , [3] a) because this 
will eventually dictate the asymptotic rate of conver- 
gence. As shown in Appendix A, we have 



because 



T{x, x) 



12r71n 



\VV{x)\\\ 



(33) 



l3Ti{x,x'](3) < E[U{x,x',(3;d) -EU{x,x',P]d)] 

< /3T2(x,x';/3), (31) 

where the functions Ti{x,x'; (3) and T2{x,x'; P) satisfy 
the relation 



1 .1 
du / dr 

"'0 



U + T — 2uT — |u — t| 



1 

12' 



T{x, x') — lim Ti(x, x'; f3) — Km T2{x, x ; (3) 



= / du [ dr — 

Jo Jo 



T — 2UT ~ \u 



xV^^^[Xriu)]V'-^^[Xr{T)]. 



Trotter composing the relation 1)31(1 n times and notic- 
(32) ing that 0(/3^) eventually contributes a term decaying 
as fast as 1/n^, it is but a simple task to establish the 
identity 



f3' 



3=0 



PA I P 

Po 2:1,0:2; 



n + 1 



Pn-j-i 2:2,2: ; 



/. {n- j)f3 



n + 1 



xTi ( 2:1,2:2; < p(a:,2:';/3) - p^^{x,x';P) < ^ .^a Y] i '^^i / '^^'i pfA {x,xi] 

\ n -\- i / z[n+L) — ; ./n .10 \ n+i 



PA I P \ PA f / (" ^ 3)P 1 m 1 
xpo 2:1,2:2; — — Pn-,-1 2:2,2: ; — ;— I T2 I ^1,2:2; 



n+1 



71+1 



J=o 

a 

n + 1 



(34) 



with the understanding that p_i (2:, 2:'; 0) = ^(2:' — x). Now, notice that in the sense of distributions, we have 

The above inequality is valid to the order of 0{(3^ /n?). p . p. , 

hnipo (2:i,2:2;/3)Ti (xi,2:2;/3) = Inn Po (2:1,2:2;/?) 

xr2 [xi,X2]l3) = 5{xi - 2:2)7(2:1,2:1). 
Multiplying it by 2(n -|- if / 13^ and using the previous 
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observation, the formula (|34|l becomes 



— - > / da;ip 1 [x,xi; — — 



< 



2(n+l)^ 

/33 



r 



in the hmit that n is large. Again in the same limit, the 
Riemann sum from the above expression transforms into 
an integral on the interval [0, 1] and combining everything 
we obtain the following theorem 

Theorem 4 

lini ^^" + ^^' [pix,x';f3)-p^^{x,x';P)] 



fc2 rl 



X 



^VVe 



12mo Jo 

It is convenient to write Eq. H35|) as 



x'^jde (35) 



p{x,x';(3)^ p^^{x,x';(3) + 

r-l 



24mo(n+ 1)2 
\e-'^"\\VVfe-^'-'^^"\ x') AO. (36) 



The d-dimensional version of Theorem 4 can be formally 
obtained by replacing ||VV^||2/mo with 



mo,i 



_d_ 

dxi 



,Xd) 



Finally, we turn our attention to the convergence of 
the RW-DPI method. It was previously proved |33| that 
for the class of potentials considered in this section the 
density matrix and the partition function of any par- 
tial averaging method is convergent to the correct result. 
However, this might not be true of the primitive and the 
reweighted methods, as well as of the standard DPI meth- 
ods. Indeed, it is known that the non-averaged methods 
suffer from what is called "classical collapse" for poten- 
tials with negative coulombic singularities [291133. |35. 36], 
for which the partial averaging method is, however, con- 
vergent. For such systems it happens that the n order 



partition functions of the primitive, reweighted, and stan- 
dard DPI methods are always -foo, yet the true quantum 
partition function is finite. This situation can be pre- 
vented by requiring that the classical partition function 
be finite. For instance, for the case of the primitive ran- 
dom series the Jensen's inequality implies 

Z^'^iP) = -jL= ^ dx J^^ dP{a) exp { - /3 V[x 

-h(7 VafcAfc(M)]dw| < -jl= f dx f dP{a) 
v27r(T^ Je Jn 

„1 n 

X / dwexpi — f3V[x + a ak^k{u)] \. 
•^0 ^ k=i ^ 



By changing the order of integration, one ends up with 



V2^ 



-l3V(x)^^ = Zel(/3) < oo, (37) 



which proves our assertion. The inequality (I37|l holds 
for the reweighted methods and the standard DPI meth- 
ods, too (for the latter techniques one uses the condition 
Wi = 1 and the discrete analog of the Jensen's in- 
equality). In this paper, the condition Zci{P) < oo is as- 
sumed to hold any time one deals with the non-averaged 
methods. 

Going back to the asymptotic convergence problem, 
we may follow the reasoning for the partial averaging 
method provided that we interpret E' to mean the aver- 
age against the Gaussian measure 



dp{z) 



^e-^'/'dz. 
/2^ 



By Jensen's inequality one proves the inequality 
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Pfp{x,x';P) 



> exp \ -P 



~ / dp{z) exp < — 13 i V[xr{u) + zTi^{u)]&u 



Pfp{x,x';P) ' 



dii{z)V[xr{u) + z ro(u)]dM > = 



Therefore, p^^ {x,x'; P) > p^^^ {x , x' ; P) . Moreover, the 
foUowing analog of Eq. H3U|) holds 

p^'^ix, x'; p) - p^^ix, x'; P) = p^^{x, x'; P) 
X I ^E' [U'{x, x', P; z) - E' U'{x, x' , P; z)f 

+^E' [U'{x, x', P; z) - E' C/'(x, x', /3; z)f + 

where we now define U' {x , x' , P; z) = V[xr{u) + zro{u)]. 
As discussed in Appendix A, we have 

PTl{x,x';P) < E'[U'{x,x',P;z)-E'U'{x,x',P;z)f 
< PT^ix,x';P), (38) 

where the functions Tl{x,x';P) and T2{x,x';P) satisfy 
the relation 



T'(x, x') = lim TUx, x'; P) = lim TL(x, x'; P) 

fc2 [\ i-l 

= — / du dr Vu(l-w)T(l-r) (39) 
"^o Jo Jo 

xV^^^[Xr{u)]V^^^[Xr{T)]. 



We also have 



T'{x,x) 



64r7io 



\VV{x) 



because 



1 ri 

du / dr \/ u(l — u)t(1 — t) = — . 
,/n ^ ^ ' ^ ^64 



We leave it for the reader to rework the previous argu- 
ments for the partial averaging case and show that for 
large n we have 

lim^^^^±^[p^^ix,x';P)~-p^^ix,x';P)] 



64too Jq 



g-e/3H||y^||2g-(i-e)/3ff 



x')dd 



Since 7r^/64 > 1/12, the previous result demonstrates 
that for n large enough p^^{x,x';P) > p{x,x';P), so 
that the convergence of the RW-DPI is eventually from 
above. Combining with Theorem 4, one obtains 



Theorem 5 

2(n+l)2 



lim 

n — *-oo 



p^ 
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[pix,x';P)-p^'^ix,x';p)] 



4mn 



x''^de. 



As for the partial averaging case, the statement of The- 
orem 5 can be written in the short form 



p{x,x';P)^pi^'^{x,x';P) 



8'mo{n + 1) 
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de. 



(40) 



From the Theorems (4) and (5) and by using cyclic in- 
variance, one easily proves the following relations: 



Corollary 1 

ZiP)-Z^^iP) 



zip) 



h^P^ J^p{x;P)\\VV{x)rdx 
24TOo(n-t- 1)2 J^p{x;P)dx 



and 



tRW 



iP) ~ Z{P) 



8mo(n-|- 1)2 1^16 ^ 3 
J^pix;p)\\VV{x)rdx 



/m P)dx 



Observation 1 We have 7rVl6 - 1/3 w 0.284 < 1/3, 
so one may be tempted to say that the reweighted tech- 
nique converges at a faster rate than the partial aver- 
aging method. However, as previously mentioned in the 
text, both the n-order LCPI and DPI reweighted tech- 
niques actually uses 2n + 1 random variables to param- 
eterize the paths. If the convention of denoting the or- 
der of an approximation by the number of variables used 
to parameterize the paths is obeyed, then the constant 
7r^/16 — 1/3 should be increased four times. In this case, 
we have 4 • 0.284 = 1.134 which means that the partial 
averaging is about 1.134/(1/3) = 3.4 times faster than 
the reweighted technique. 

Observation 2 The asymptotic relative errors for the 
partition functions shown in Corollary 1 can be evalu- 
ated during the Monte Carlo procedure if so desired. It 
is a fact established in several occasions l.'i. 'AO that 
the convergence of the partial averaging density matrix 
and partition functions for all series representations is 
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monotonically from below. In particular, the PA-DPI 
subsequence n = 2^ — 1 has the same property since it is 
identical to the respective subsequence of the PA-LCPI 
method. However, it might be possible that the partition 
function for the reweighted methods are monotonically 
decreasing from above for the n — 2^^ — 1 subsequence. 
In fact. Golden 'sT] and Thompson 38] have shown that 
the partition function for the subsequence n = 2'^' — 1 of 
the trapezoidal Trotter DPI method is monotonically de- 
creasing and this might be true of the RW-DPI method, 
too. 

Let us remember that there are potentials, as for in- 
stance the potentials with negative coulombic singular- 
ities, for which the non-averaged methods do not con- 
verge. Conversely, there are smooth and bounded from 
below potentials, as for instance V{x) = exp(a;^) for 
which the non-averaged methods are convergent to the 
correct result yet the partial averaging method is not con- 
vergent because V{x) — exp(x^) does not have a finite 
Gaussian transform. For such potentials, it is expected 
that the Theorem 5 as well as the second part of the 
Corollary 1 are still true. 

We shall reinforce the conclusions of this section by 
verifying the theorems (4) and (5) for the simple case 
of the quadratic potential V{x) — mouo'^x'^ /2. Again 
we use atomic units and set mo — \^ u) — 1, and 
/3 = 10. The evaluation of the n-order partial averag- 
ing and reweighted elements p^^(0;/3) and P^{0;(3) is 
analyzed in Appendix B. As discussed in Ref.'Ts, for each 
method Alt the convergence constant 

CMt = lim — —2 

can be obtained numerically by analyzing the asymptotic 
slope of the equation 

= (4n + 2)2(n+ 1/2) K0;/3) -pr+2(0;/3)] , 

as a function of n. More precisely, we have CMt = 
lim„^oo c^^* — Cn-i- the other hand, with the help of 
the exact density matrix p{x,x';P) of the quadratic po- 
tential [s^ and of the relations H36|) and H4U|) . one com- 
putes 

^4 JO 

xp[x,0; (1 - 61)^] x"^ = 0.0713 

and crw — — [(37r^/16) — 1]cpa = —0.0606, respectively. 
The plots in Fig. 01 show that indeed, the current slopes 
c^^*— c*{^]^ converge to the corresponding values predicted 
by the theorems (4) and (5). 

To conclude this section, we analyze how smooth real- 
istic three dimensional potentials must be to fit the hy- 
pothesis of the the theorems (4) and (5). A prototypical 
example is the tridimensional spherical potential 

= + o<a<i, (41) 




FIG. 4: The current slopes c*^' —Cnli (solid lines) are shown 
to converge to the values predicted by the theorems (4) and 
(5) (dotted lines). 

for which both the partial averaging and the reweighted 
DPI methods are convergent because the potential V{r) 
is a Kato-class potential having a finite Gaussian trans- 
form and because Zci{f3) < oo. The reader may easily 
verify that ||Vl/r"|p = a^/r^°'+'^ is locally integrable if 
and only if a < 1/2. Therefore, if a < 1/2, the theorems 
(4) and (5) apply and the convergence of both methods 
isO(l/n2). 

On the other hand, if a > 1/2, the convergence cannot 
be 0{l/n^) because the convergence constants are -l-oo. 
This can be proved by using the additional information 
that the density matrix for the Kato-class potentials is 
continuous and strictly positive. In particular, there is 
e > and > such that p{x,P) > e for all r < 77 < 
1. Therefore, looking at the bounds for the partition 
functions given by Corollary 1, we have 

[ p{x;(3)\\\'V{x)fdx>ATTea^ [\-^°'dr 

> Anea^ / r~^dr — +00. 
Jo 

We have treated this problem explicitly in order to show 
that the nature of the singularities of the potential af- 
fects the rate of convergence even if the singularities are 
oriented "upward." Therefore, in "pushing" the Monte 
Carlo simulation to the limits, the reader may want to 
actually remove these singularities if they are physically 
irrelevant. He/she can do this either by a simple trun- 
cation or by approximating the singularity with a better 
behaved one. 

Other prototypical examples of potentials are those 
having negative singularities 

0<«<1. (42) 

For such potentials, the classical partition function is 
not finite and the reweighted technique does not prop- 
erly converge. However, the partial averaging method is 
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convergent and, if a < 1/2, the asymptotic convergence 
is 0(l/n^). The findings of this section demonstrate that 
"smooth enough" potentials may actually be discontinu- 
ous in the three dimensional space. 

IV. SUMMARY AND DISCUSSION 

A central theme of the present paper has been the char- 
acterization of various path integral approaches and the 
exploration of their interconnections. The Kag interpre- 
tation of the Feynman approach is a valuable tool for 
such an analysis. We notice that it is difficult and un- 
natural to introduce the random series representation by 
means of the Trotter product rule. Indeed, in order to 
show that the path integrals 

/ V[xr{u) + Bl{a)]du (43) 
Jo 

are correctly defined, one utilizes the fact that, with prob- 
ability one, the Brownian paths are continuous. This 
property of the Brownian motion is not readily avail- 
able from the Trotter product rule. However, as we have 
shown in Section II, the discrete methods can be directly 
derived from the Feynman-Kag formula by simply re- 
placing the integrals given by Eq. I|43|) with appropriate 
quadrature sums. 

We have explored at some length two particular imple- 
mentations of path integral methods: the Levy-Ciesielski 
approach and the associated DPI technique. We have 
considered primitive, partial averaged and reweighted 
versions of this methods. As discussed in Section III, 
the Levy-Ciesielski approach is of particular importance 
because its n = 2*^ — 1 subsequence can be rationalized 
both as a series and as a discrete method. This dual 
character is valuable for several reasons. For example, 
it provides a convenient and rigorous reformulation of 
Coalson's findings linking series and discrete path inte- 
gral methods, and, as illustrated by Eq. lf77|l . suggests 
a means for reducing the numerical overhead associated 
with path construction. Using the unified framework 
the Levy-Ciesielski approach provides, we have shown 
that the methods introduced by Kole and De Raedt 
for systems with negative coulombic singularities as well 
as those introduced by Titantah, Pierleoni and Ciuchi 
[3^ for the polaron problem are discrete versions of the 
partial averaging approach. Furthermore, Theorem 2 of 
Section III suggests that these previous methods can be 
implemented in a more robust manner using the Levy- 
Ciesielski series approach. 

We have been able to characterize the convergence 
properties of the partial averaging and reweighted DPI 



approaches and, therefore, of the n = 2^ — \ subsequence 
of the corresponding LCPI techniques. In this respect. 
Theorems 4 and 5 of Section IV provide sharp estimates 
of the convergence constants for the calculation of den- 
sity matrix elements for both the partial averaged and 
reweighted methods. To our knowledge, this is the first 
time that such exact convergence constants have been 
established. Beyond their intrinsic interest, knowledge 
of these convergence constants can be used to devise an 
improved numerical implementation of the Feynman-Kag 
approach. In particular, the results of Section IV indi- 
cate that the convergence constants for the reweighted 
and partial averaged methods are related by the formula 
crw — ^[(37r^/16) — V\cpA for all pairs of points (x,x') 
and for all /3 > 0. Because the leading terms in l/v? thus 
cancel, the approach defined by the equation 

, , ^, pT{x, x'- /3) + [(3^V16) - l]plHx, x'- (5) 



has an asymptotic convergence better than ©(l/n^), i.e. 
hm -^^^^^ [p(x, x'- /3) - pUx, x'; /?)] = 0. 

In fact, we believe that if the potential V[x^ has also a 
well behaved second derivative, the convergence order of 
the new method is ©(l/n^). 

Finally, we note that with the help of Theorems 4 and 
5, the asymptotic behavior of the so-called T-method and 
H-method energy estimators (c.f. Section IV of Ref. IT^ 
can be examined. In particular, it should be possible 
to deduce the convergence constants for these estimators 
from those of the corresponding density matrix expres- 
sions. We leave a detailed analysis of such issues for 
future discussion. 
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APPENDIX A 

It is well known that A,B > 0, and a = C jyj AB 
such that |q;| < 1, then the following Mehler's formula 
^b] holds for all / and g whose squares have finite Gaus- 
sian transforms 
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rrr / n f , f , ^ 1 f I B + y^A ^ 2xyC\ ^, 

= fdxfdy^^=L=e^p(~l-^^:f^p^)f{xo + xVA)g{yo + yVB) (Al) 
Jr 7r 27r ^1 - a2 V 2 1 - y 

= ^E"' /dx/ d2/e-(^'+«')/2^fe(x)i/fe(y)/(xo + a;Vl)g(yo + 2/VB). 



In the above, the functions (x) are the normahzed Her- 
mite polynomials corresponding to the Gaussian weight 



d/i(x) 



1 



2tt 



They form a complete orthonormal basis in the Hilbert 
space i^(M), which is endowed with the scalar product 



(^10) = / ^P{x)(t>{x)d^iix). 
Js. 

Let us notice that according to our hypothesis, the func- 
tions f{xo+x^/A) and g{yo + x\/B) as functions of x are 
square integrable against dfi{x) and thus they lie in the 
Hilbert space ^^(R). 

By repeated integration by parts, the formula I|A1|I is 
shown to equal 

°° C''-(fc) 



where in general fj^ (xo) is the fc-order derivative of 

7Axo)= I -i=e-^^/(2-4)/(xo + ^)dz. 
JR v27r/l 

Let us notice that the series ljAl|l can be extended to the 
case (3 = 1, too. Indeed, the last series in Eq. (|A1() for 
the case a = 1 is nothing else but the Bessel series 



(Hk\f{xo + -VI)) {Hk\g{yo + -Vs; 

k=Q 

which is convergent to 

(/(xo + -VI)|.9(yo + -\/B)) 

= / f{xo + x\/A)g{yo + xy/B)dfi{x). 



Next, we proceed to establish the inequality (|31|l from 
section III.C. We start with the identity 



E [U{x, x',P;a)-E U{x, x' , /3; a)f 
E U{x, x', P; a)^ - [E U{x, x' , /3; a)f 



(A3) 



Clearly, we have 

E U{x, x', (3; a) = VuAxr{u)]. (A4) 

Moreover, 

EU{x,x\P] a)^ 

^ [ du [ dTEV[Xr{u) + (TB°]V[Xr{T) + aB°] (A5) 

and the variables B^ and B^ have a joint Gaussian dis- 
tribution of covariances 

EiB',f = u{l-u), E(SO)2=r(l-T) 



and EiB'^B'l) 



u + T — 2uT — |u — r| 



This covariance matrix is independent of any particular 
representation of the Brownian bridge and therefore can 
be computed with the help of any basis. For instance, 
using the Wiener- Fourier basis, the last term of the above 
formula reads 

^KB'r) = ^f: (A6) 



fc=i 



and the sum of the above series can be shown to equal 
{u+t~2ut—\u—t\)/2. It is useful to define the quantities 

Go{u,t) = (T^E{B^^B°) and 



A'o{u,T) = rl{u)Tl{T)-Go{u,Tr 



Then, 



EU{x,x',P;dy = / du dr dx dy 



X exp 



10 Jo JR JR 27rAo(M,T) 
lx^Tl{T)+y^Tl{u)~2xyGo{u,T)- 



2 Aliu,r) 

xV[Xr{u) + x]V[Xr{T) + y]. 

Using the expansion ljA2p . one may write the above 
integral as the sum of the series 

EU(x,x',f3;af ^ Y.T\ / dTGo(u,T)'= 



X vi%Au)]V^rl[Mr)], 
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if^\ 

where q{x) is the k order derivative of Vu,o{x). With 

the help of Eq. ljA4|l . one recognizes the first term of 

the above series to be [KU{x,x' ,P; a)]^, so that Eq. (|A3|) 

becomes 



E 



[U{x, x\ /3; a) - E U{x, x',13; a)f = ^t 



k=l 



kl 



X l\u[\TGo{u,rfv[%r{u)]V^rl[xr{r)]. (A7) 



Now, we make an important observation: as its eigen- 
function expansion ljA6p shows, Go{u, r) is a positive def- 
inite integral kernel L^([0, 1]) i^([0, 1]) and it is not 
difficult to verify that all Go{u,t)'' are positive definite. 
Therefore, 



'd.. rdTGo{u,T)''vi%r{u)]vi%r{T)] > 



for all fc > 1. Considering only the first term in the series 
(|A7|) . we obtain the inequality 

E[U{x,x',(3;a) - EU{x,x',l3;a)f > f3Ti{x, x' ; (3) (AS) 



where 



Ti{x,x';P) 



s=2 rl 



du dTE[B°B°] 
Jo 



mo Jq 



It is not difficult to see that as /? ^ we have ^o{u) — > 
and so, 



T(x, x') = lim Tiix, x'; (3) — 

/3^o mo 

X f\u f\TE[B°B",]v(^^\xr{u)]V^^\xr{T)]. (A9) 



To prove the second inequality in Eq. (|31|l . one uses 
the inequality 1/kl < l/(fc — 1)! and the positivity of the 
terms of the series l|A7|l to establish the inequality 



OO „i „i 

E[U{x,x\P;a)^EUix,x',P;a)f <Y,-n / du dTGo{u,T)''+'vi'o'\Mu)Wi'o'\Mr)] 

k=0 ■ "^0 



fe=0 



1 



27rAo(u,T) 1^ 2 A5(u,t) J 

I 



Therefore, 

E[U{x,x' , I3]a) - EU{x,x' , (3]a)f < PT2{x,x'; (3), 
where 

T2{x,x']P) = — [ du [ dTE[BlB°] [ dx [ dy 
mo Jo Jo Jm Jr 

lx^rl{T)+y^rl{u)-2xyGo{u,T) 



The relations (|38|l and H39|) follow by a similar reasoning 
and their proof is left to the reader. We only mention 
that one starts with the fact that the series HA1|I is well 

(AlO) 

defined and convergent for a = 1 too, as shown in the 
beginning of the present appendix. 



27rAo(M,T) [ 2 



A§(u,r) 



APPENDIX B 



xV''^'^[Xriu) + x]V^^^[Xr{T) + y] 

Again, as — + we have Fq (m) and 



T(x,x') = lim r2(x, a;'; /3) 



In this section we discuss the computation of the 
matrix element (0|e^'^-^|0) for the quadratic potential 
V{x) = mouP'x^ 12 by means of the standard DPI method 
and of the partial averaging and the reweighted DPI 
^1 /•! methods. The density matrix for the quadratic potential 

X du drE [B°B°]V'-^^[xr{u)]V'^^^[xriT)]. (All) is known analytically (see Ref.'s^ and we do not repro- 
° ° duce it here. For a standard DPI method specified by the 

The relations (|A8(l . ljA9|l . (|A10(l . and ljAll|l combined quadrature points = uo < ""i < • ■ ■ < < "n+i — 1, 
prove the equations (|31|l and H32|) from Section III.C. by the increments 9i — Ui+i — Ui, and by the weights 
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Wo,Wi, . . . ,Wn+ 



1, the formula (|12|l becomes 



where the tridiagonal matrix A is defined by the relations 



'XPa^e2{^2 -Xi).. .p„2g^{Xn - Xn-l)p„2g^{Xn) (Bl) 

i=l 



X exp ■ 



Remember that = l3fi^/mQ. If we set x'^ = 
{xi,X2, • ■ • , Xn), the above n-dimensional integral can be 
written in the compact form 



1/2 



11 Oirrr 



1/2 



\i=0 



det 



g-2-As/2j- 



27r 



(B2) 



where the matrix A is the tridiagonal matrix defined by 

A^ 

and 



1 1 2 
\ :t— + mQU) f3wi for 1 < i < n 



^i,i+l — Ai^i^i — 



for 1 < i < n — 1 . 



The values of the quadrature points and the correspond- 
ing weights for the trapezoidal rule are well known, while 
for the Gauss-Legendre quadrature scheme the reader 
may use the routine given in Ref. l|2nt ). 

The zero order partial averaging density matrix for the 
quadratic potential has the explicit expression 



PO {x,x']P) = Pa-^ix - x) 



X exp 



6 



(B3) 



Using Eq. (|B3|I . the reader may easily deduce that the 
corresponding rt-order PA-DPI density matrix is 



pr(0;/3) 



X exp 



f n+l 
12(n+ 1) 



(n+l)/2 



det 



1-1/2 



n + 1 moUj'^P 



3(n+l) 



for 1 < i < 



■ ^ Wix'fdu I . and 



Ai^i+i — Ai^i,i 



n + 1 niQbj'^jB 



cr2 6(n + l) 



for 1 < z < n — 1. 



Finally, the zero order reweighted density matrix has 
the form 



pI^{x, x'; 13) = {x' -x){l+ 



X exp 



9 /2 / 

X + X + XX 



(B5) 



TT^ l3^h^L0^{x + x') 



l\2 



W 1 + l3'^h?Lu^/6 



which can be deduced by direct integration. Let us set 



2 /32ft2^2 



6(n+ 1)2 



Then, 



pr(0;/3) 



1 1 



\ 277 tr^ 7^2 



-1/2 



(B6) 



where the tridiagonal matrix A is defined by the relations 



Ai i — 2 



n + 1 ^ mQu'^P I" \2 TOo^^/3^w* 



1) V 



for 1 < i < n and 



Ai i+i — Aij^i i — — 



3(n + l) Vl6/ 77,2(^ + 1)3 



n+1 ( ^ \ rrifih^ fi'^ lo'^ 

cr2 ^6(n+l)~vT6/ r;2(n + 1)3 



, (B4) for 1 < i < n - 1. 
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